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We perform quantum Monte Carlo calculations to determine how the Renyi entropies, S n , of the 
interacting Fermi liquid depend on Renyi order, n, and scale as a function of system size, L. Using 
the swap operator and an accurate Slater- Jastrow wave function, we compute Renyi entropies for 
spinless fermions interacting via the Coulomb and modified Poschl- Teller potentials across a range of 
correlation strengths. Our results show that interactions increase the Renyi entropies and increase 
the prefactor of their scaling laws. The relationships between Renyi entropies of different order n 
are also modified. Additionally, we investigate the effect of the swap operator on the Fermi liquid 
wave function to determine the source of the L log L scaling form. 



Entanglement is a central concept in quantum infor- 
mation theory, and has become increasingly important in 
studies of quantum many body models of condensed mat- 
ter pJ|Z] • The Renyi entropies computed from a spatially 
reduced density matrix provide a measure of entangle- 
ment, through non-local correlations in the ground state 
wave function [2j |3j [8]. The scaling laws for these en- 
tropies provide a "fingerprint" for quantum phases [6] [9]- 
[T2] . For ID systems these scaling laws can often be com- 
puted analytically using results from conformal field the- 
ory [12] [13] . In higher dimensions the Renyi entropies are 
more difficult to compute and are often calculated using 
numerical techniques [T4H25] . No general results for in- 
teracting systems in higher dimension are available, mod- 
els must be addressed individually and results can differ 
greatly from their non- interacting counterparts [THJ [19] . 

One of the most fundamental models in condensed 
matter is the Fermi liquid [26] . The Renyi entropies of the 
non-interacting Fermi liquid in two dimensions have been 
conjectured to scale as LlogL [2T j [27H3Q] . This scaling 
law is equivalent to the Widom conjecture and has been 
verified numerically for lattice and continuum Hamiltoni- 
ans [29H33] . High dimensional bosonization and confor- 
mal field theory studies have predicted that the scaling 
law of the Renyi entropies for the interacting and non- 
interacting Fermi liquid are identical [34-36 ]. These re- 
sults imply that the relationship between S n of different 
order, S n = ^(1 + should also hold when interac- 

tions are included [T3j |34] . 

In this work, we use variational quantum Monte Carlo 
(VMC) to test these predictions. Specifically, we com- 
pute the Renyi entropy scaling forms for two differ- 
ent inter-particle potentials and compare them to the 
non-interacting results. VMC provides an accurate 
framework for the calculation of interacting and non- 
interacting quantum states in a many body framework, 
and is expected to produce reliable results for well opti- 
mized trial wave functions [37]. High quality trial wave 
function forms and robust optimization schemes have 
been developed to accurately compute observables [37]- 
[40] , These have been used to calculate benchmark en- 
ergies for the electron gas and Fermi liquid parameters 



for the Coulomb potential, Helium-3, and several related 
Fermi liquid systems [37] [38 | I4T1443] . 

Hamiltonians and Trial Wave Functions- We choose 
two inter-particle potentials, the Coulomb and the mod- 
ified Poschl-Teller potential (MPT), to test if our re- 
sults are sensitive to the range of the interaction poten- 
tial. Spin polarized electrons interact via the long range 
Coulomb potential [44]. The MPT potential is a well 
tested short range potential for the Fermi gas [45-47 . If 
the range of the potential is important to the Renyi en- 
tropies then we could expect some qualitative difference 
between the scaling laws computed for these potentials. 

The Hamiltonian for the spinless non-interacting elec- 
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tron gas in Hartree units is Hfg = — ~2 • ^he 
exact fermionic ground state for this Hamiltonian is 
^fg(R) = (R\(Uk<k f c ll°))- For the interacting elec- 

V 2 —1 

tron gas H H eg = ~ J2i ~t + ^i<j r ij + c ( r s)> where 
r ij ~ \ r i ~ r j\ 1S the distance between electron i and 
j, and the constant, C(r s ), is due to a uniform positive 
background charge. For the MPT potential Hmpt — 
~ 12i + Hi<j cosh -2 (r^). The parameter V can 
be tuned to increase the correlation of the system and the 
effective range is set to the Wigner-Seitz radius, r s = 1 

051- 

We use a Slater- Jastrow trial wave function, 
J(R)^ fg(R), f° r these interacting Hamiltonians. The 
Jastrow is of the form J(R) = exp (X^<j ^( r ij)) an d 
adds additional correlation between particle pairs. It 
is a reasonably accurate form that can be quickly com- 
puted [33 HH1 [49]. For the electron gas wave function, 
^ heg(R), we eliminate the divergence in the electron- 
electron energy as r^ —> by including a cusp in the 
Jastrow [S3 EQ- For the MPT wave function, ^mpt(R), 
there is no singularity in the potential and therefore no 
cusp in the Jastrow. 

Renyi Entropy and the Swap Operator- The Renyi en- 
tropies are defined as 

S n {pA) = Y^logOW)) (1) 

with the spatially reduced density matrix, p^, of a system 
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partitioned into two regions, A and B. Using VMC we 
are able to compute unbiased estimates of S n for trial 
wave functions using the swap operator [52] . 



exp((l-n)S n ( M )) = Tr(^) 

^(g)^--. <8>#| SWAP n |*( 



(2) 



The swap operator computes Tr(//|) by cyclically swap- 
ping coordinates between n statistically independent 
copies of the system, each sampling \I/ 2 . For instance, 
we can compute S 2 using, 



<8>tf| SWAP 2 |^(8 ^ 



(3) 



where R\ represents the electronic spin and spatial co- 
ordinates in region A for copy z, R % = (R % Al R B ) an d 

Computing Renyi Entropies- In practice, the expecta- 
tion value of Tr(p A ) decays rapidly as region A increases 
in size [21]. To compute Renyi entropies for large re- 
gions we factorize the swap operator into two parts: sign, 
S a , and magnitude, Sm [2T] - These two pieces are com- 
puted separately then summed to obtain the full value. 
This factorization is not an approximation and does not 
change the values we compute for S n . The magnitude 
component can be further factorized into Sm = SV + Sjv? 
with Sn due to swapped copies having different numbers 
of particles in region 5, and Sy from the absolute value 
of the swap when they do have the same number. The 
estimators for the full factorization, 

^ log (( SWAP„)) = (S N ) + (S v ) + (S a ) , (4) 

can be found in the supplementary material. 

Source of the Area Law Violation- In Fig. [I] we plot 
the contributions to 62 from the factorization in eqn. [4] 
for the non-interacting Fermi liquid and the interacting 
electron gas at r s = 5. The contribution to S2 from the 
sign of the swap operator grows as L increases, while the 
other contributions plateau or decrease. This indicates 
that as L gets large the sign of the swap operator is the 
only component that can be responsible for the leading 
scaling behavior. We note that previous studies found 
the same qualitative behavior [2T] , 

We investigate the spatial structure of the amplitude 
of the swap operator in Fig. [2] The upper plot shows 
SWAP 2 (n) = ^(R 1 )^(R\,R 2 B )^(R 2 )^(R%R 1 B ), with 
7*1 a coordinate in Ri: R\ = {7*1, r 2 , . . . , rjv e }- This plot 
is produced by swapping all particles in region B then 
scanning a single particle initially located at n , through 
region A and plotting the amplitude of the swap operator. 
The lower plot shows the normalized average amplitude 
of the swap operator as a function of distance from the 
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FIG. 1. S2/L as a function of logL, L y/ir(N), for (Left) 
the non-interacting electron gas and (Right) the interacting 
gas at r s = 5. The bold black line is 6*2- Sn, Sa, and Sv are 
the factorization of & as in eqn. [4] The log L boundary law 
violation comes from the sign of the swap, it is the only term 
that continues to increase as L gets large. For the interacting 
system at r s = 5, the Jastrow increases the contribution from 
Sv- 
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FIG. 2. (Top) We scan a single particle from n through 
the circular region A and plot SWAP2(ri), defined in the 
Letter body, for The white dot is the scanned particle 

at it's original position, n, in the lower left of Region A. 
The swap operator exchanges the red (blue) particles from 
copy R l (R j ) outside of region A. The red (blue) regions are 
positive (negative), with saturation representing magnitude. 
(Bottom) The average amplitude of SWAP2(r), where r is the 
distance from the center of region A. SWAP2(r) is larger than 
it's average value in the interior and decays rapidly near the 
boundary. The sign behaves similarly. 



center of region A. This average amplitude, just as the 
average sign, is approximately constant in the center and 
decays as particles are moved out to the edge. This sug- 
gests that the change in sign of the swap operator due 
to particles near the boundary drives the log violation of 
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the area law. 

Heuristically, we can justify the Jastrow modifying the 
L log L and L scaling terms. The Jastrow can change the 
leading scaling law of the Renyi entropies by changing 
the particle distribution. It can also create additional L 
dependence by destroying the short range correlations in 
the Jastrow for the particles in region A near region B. 
We expect the difference in the Renyi entropies between 
the interacting and non-interacting systems to be well 
fit by these corrections, AS 2 /L = (S^ EG - S 2 FG )/L = 
I + mlogL. 

Computational Details- All QMC calculations are per- 
formed in QMCPACK[53j. We parameterize the Jastrow 
by a flexible B- Spline function with fixed cusp conditions 
depending on the potential and optimize them using a 
variant of the linear method [40] [54] . A more accurate 
Slater- Jastrow-Backflow wave function was also tested 
[55] . Using this wave function Renyi entropy results for 
6*2 at r s = 5 were the same as the Slater- Jastrow values. 
The only changes were in the orbital contributions to the 
Renyi entropy. Backflow increased the determinant and 
decreased the Jastrow contributions to the swap magni- 
tude. 

For the electron gas, the swap operator for S2 was used 
without factorization for circular regions A up to (N) = 
16 for a simulation cell of 137 electrons and with the sign 
factorization technique for (N) = 16, 25, and 36 in 261 
electron systems. We computed 62 for regions of (N) = 
1 though 16 for the 137 particle simulation cell using 
the unfactorized estimator for the MPT potential. When 
factorization is unnecessary we are able to increase our 
computational efficiency by computing the swap operator 
for several regions sizes at the same time. Approximately 
10 6 CPU hours were spent on the smaller systems and the 
same on the larger. For the higher order Renyi entropies 
we were restricted to (N) = 16 or fewer electrons for the 
261 electron system. 

We rescale all lengths, L — ^ L/(^r s ) = ^7r(N), so 
that S n (L) for the non-interacting Fermi liquid has no 
r s dependence. Under this rescaling, all interacting S n 
will collapse to the non-interacting S n if interactions are 
irrelevant. We use a logL x-axis to highlight the log 
linear nature of S n /L and AS n /L. 

Results- We find that interactions increase the Renyi 
entropies and increase the prefactor in their scaling law. 
More strongly correlated Hamiltonians produce ground 
state wave functions with stronger particle correlations 
and larger Renyi entropies. We also find that Coulomb 
interactions change the simple relationship between the 
Renyi entropies of different order n. 

The non-interacting data is well fit by the leading scal- 
ing form predicted by the Widom conjecture, S(L) = 
m (L/l ) log(L// ) with ra = 0.032(2) and l = 0.113(5) 
[13 [301133]. 

As shown in Fig. [3] Coulomb interactions modify 
the Renyi entropies when the Coulomb energy becomes 



r s 


m 


I 


Z 


Corr 


1 


0.005(1) 


0.001(1) 


0.922(1) 


0.991 


5 


0.035(1) 


0.001(3) 


0.614(1) 


0.990 


10 


0.053(2) 


0.027(3) 


0.419(1) 


0.988 


20 


0.073(2) 


0.065(4) 


0.272(1) 


0.985 



TABLE I. r s is the Wigner-Seitz radius of the electron gas. 
r s — is the non- interacting system. AS2/L = l + m\ogL. Z 
is the quasi particle renormalization factor. Corr= (Evmc — 
Ehf) I \Edmc — Ehf) is the ratio of correlation energy the 
wave function recovers and provides a metric for wave function 
quality. 
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m 


I 


Z 


Corr 


10 
20 


0.032(1) 
0.062(1) 


0.007(2) 
0.028(1) 


0.708(3) 
0.487(4) 


0.984 
0.981 



TABLE II. Vb is the strength of the potential for the modi- 
fied Poschl- Teller (MPT) potential. All other labels are as in 
Table H] 

larger in magnitude that the kinetic energy. At r s = 1 
the Coulomb energy is approximately equal to the ki- 
netic energy and interactions do not significantly alter 
S2 • As Coulomb interactions become stronger, the Renyi 
entropy 62 increases. For r s = 20 the Coulomb energy 
is an order of magnitude larger than the kinetic energy, 
and the Renyi entropy is increased relative to the non- 
interacting Fermi gas. The higher order Renyi entropies, 
reported in the supplementary material and shown in Fig. 
[4] show the same qualitative trend. 

We find AS2/L is well fit to the two parameter form 
I + m log L. The fitting parameters for the Coulomb po- 
tential at all r s are presented in Table [I] and for the MPT 
at all Vb in Table [IT] The Coulomb fit is plotted in the in- 
set of Fig. [3] As the strength of the inter-particle poten- 
tial is increased so too is the prefactor of the leading scal- 
ing term, m. We compute the quasi particle renormal- 
ization factor, Z, to determine how strong inter-particle 
correlations are and compare scaling laws between differ- 
ent model potentials [49]. The scaling laws for the MPT 
potentials show qualitative agreement with the interact- 
ing electron gas with the same Z. As Z increases so does 
the slope of AS2/L. These slopes are almost linear with 
Z and are plotted in Fig. [5] We note that finite size 
corrections to Z may not preserve this relationship [49 . 

On the left of Fig. [IJ AS n /L is plotted for the Coulomb 
potential for a region of size (N) = 6.25. The simple rela- 
tionship between non-interacting S n does not hold when 
Coulomb interactions are included. Interactions increase 
the Renyi entropies for all n. As shown on the right of 
Fig. jij upon rescaling, ^ — > f^^f, the non interacting 
Renyi entropy is made constant. The difference between 
the interacting and non-interacting S n decreases as n gets 
larger. The growth of the error bars as L and n increase 
prevent a good fit for their scaling laws and the relation- 
ship between S n of different order n. 



4 




0.6 1.0 1.4 1.8 2.2 

logL 

FIG. 3. (Main Plot) S2/L for the spin polarized interacting 
electron gas at r s = 1, 5, 10, 20 and the non-interacting Fermi 
liquid. L is scaled so that without interactions all lines lie 
on top of f^G results. The high density r s = 1 data falls on 
top of the non-interacting liquid while the lower density, more 
strongly correlated, gasses have larger Renyi entropies. The 
^ fg reference is the exact non- interacting S2 computed using 
the same number of particles as the interacting system for 
each L. (Inset) AS 2 /L = (S* EG - Sf G )/L plotted against 
logL. The Renyi entropy for the interacting case appears 
to scale the same as the non-interacting case with a larger 
prefactor. The lines are from a fit to AS2/L = I + ralogL 
shown in Table HJ 
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FIG. 4. For the Coulomb potential: (Left) AS n /L, for n = 
2,3,4 for region A size logL « 1.5, (N) = 6.25. S n increases 
as correlations increase and decreases as n increases. (Right) 
t~t— ^r- for n = 2, 3, 4. Besides some small L oscillations, this 
scaling makes the non-interacting Renyi entropies linear in n. 
Interactions modify this relationship and the rescaled Renyi 
entropies decrease as n gets larger. 



Error analysis- For these calculations there are two pri- 
mary sources of error, finite size effects, and trial wave 
function bias. Finite size effects come from two places, 
one is a function of the number of electrons in the system, 
the other is from subleading corrections to the Renyi En- 
tropy scaling law. By integrating out regions that are a 
small fraction of the system size, we minimize interac- 
tions of region A and it's periodic neighbors. For small 
subregions we see oscillations, small relative to S2, that 




1 - Z 

FIG. 5. Slope difference, m (AS2/L = I + ralogL), as a 
function of 1 — Z, one minus the quasi particle renormalization 
factor, for the Coulomb and modified Poschl- Teller (MPT) 
potential at all correlation strengths considered. Lines are 
guides to the eye from linear fits constrained to go through 
the origin. 



decay as the subregion length increases for both the inter- 
acting and non- interacting cases [33] . Fits are performed 
on AS2/L to minimize the effect of oscillations on the 
fitted parameters. 

The degree to which the trial wave function represents 
the Fermi liquid depends on it's overlap with the exact 
ground state wave function. We estimate the quality of 
our trial wave function by looking at the amount of cor- 
relation energy the VMC recovers relative to the fixed 
node diffusion quantum Monte Carlo result, E corr — 
Edmc — Ehf [37]. We recover more than 98% of this 
correlation energy for all densities and potentials, shown 
in Table [T] and [Til Also because the Slater- Jastrow and 
Slater- Jastrow-Backflow Renyi entropies are similar, we 
expect them to be accurate. 

Conclusions and Future Work- In this work we have 
performed the first ab initio calculation of the Renyi en- 
tropy of the interacting Fermi liquid. Our results show 
that, for this system, interactions modify the scaling be- 
havior of the Renyi entropies. As particle correlations, 
quantified by Z, are strengthened, the prefactor for the 
leading scaling behavior of the Renyi entropies increases. 
These results hold for all wave function forms, correla- 
tion strengths, and inter particle potentials considered in 
this work. We also find the relationship between S n of 
different order n changed. The source of the disagree- 
ment with recent theoretical predictions is not clear and 
deserves further investigation [34-36J. Work is currently 
underway improving the efficiency of computing higher 
order Renyi entropies and investigating related systems 
such as the paramagnetic electron gas and the Wigner 
crystal. 
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SUPPLEMENTARY MATERIAL 



VMC Estimator 



The VMC estimator can be factorized as, 

^2 = -l0g((* 2 



^ 2 

log ((S# 2 J) 



SWAP 2 



* 2 



log (Sn u n 2 ) x 



|SWAP 2 | 

(^iVi,JV 2 ) 



(S N ) = 
(S V ) = 



(S N ) + (S a ) + (S v ) ■ 
-\og((5 NltN2 )) 

|SWAP 2 i 




(s# 2 ; 
swApal 



(5) 
(6) 

(7) 

(8) 
(9) 

(10) 
(11) 



where Sn is the number fluctuation contribution, S a is the sign contribution, and Sy is the magnitude contribution. 
The Kronecker delta, 5jvi,jv 2 ? is one when copy 1 and 2 have the same number of particles in region B and zero 
otherwise. 

To more efficiently compute the sign we sample a different probability distribution, the swap magnitude, 



(Sa) 



V(R 2 ) 



^(R 1 ) V(R 2 ) 



(12) 



where R\ represents the electronic spin and spatial coordinates in region A for copy z, R 1 = (R % A ,R % B ) an d 
^(R) = (R\^)- This computes the sign of the swap operator from the regions in R2 that contribute most to 62- 
The magnitude and particle number contributions decay more slowly than the sign contribution and are sampled 
over the original distribution, \I/ 2 \I/ 2 . 



Higher Order Renyi 



Because the expectation value of the swap operator can be negative when the value is small and the error is large, we 
are unable to plot the full range of values for £3 and S4 in Figure [6] It can be seen that for the region sizes we are able 
to compute, the qualitative trends for S n (L) shown for £2 in the Letter also hold for £3 and 64, S^ EG (L) > S EG (L). 



